Relaxation and Rheology in Dense Athermal Suspensions 
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We study relaxation and rheology of dense athermal suspensions of frictionless particles close below 
the jamming density. Our key quantity, the relaxation time — determined from the exponential decay 
of the energy after the shearing has suddenly been switched off — is argued to be a determining factor 
behind the algebraic divergence of various quantities as the jamming density is approached from 
below. We also define and measure the "dissipation time", which is obtained directly in shearing 
simulations and find that it behaves similarly to the relaxation time. Comparing shear viscosity 
with the expression for the dissipation time we identify a non-divergent factor that explains the 
need for correction terms in the scaling analyses of the shear viscosity. 

PACS numbers: 63.50.Lm, 45.70.-n 83.10.Rs 



As the volume fraction increases in zero-temperature 
collections of spherical particles with repulsive contact 
interaction, there is a transition from a liquid to an amor- 
phous solid state — the jamming transition. This transi- 
tion has for quite some time been studied through simula- 
tions in two different ways: by examining static packings 
generated by compressing and relaxing random packings, 
and by driving the system with a shear deformation. 
Whereas it was first commonly expected that these two 
approaches would show the same behavior, the evidence 
now suggest that they are clearly different. One example 
is the difference in the behavior of the pressure above the 
jamming density, 0j, ~ (</) — which is linear, 

?/ = 1 for static packings [1] but appears to be y « 1.1 for 
the shear-driven case [2]. Another example is the isolated 
mode in the spectrum that dominates the behavior of the 
shear-driven system close to the transition [3] but which 
is not present in static packings. 

One way to study the shear-driven transition is to try 
and eliminate the complications related to the softness 
of the particles and instead try and determine the be- 
havior of hard particles. This is usually done by driving 
with sufficiently low shear rates, 7, such that the parti- 
cle overlaps become negligable — this is the linear region 
where many quantities are linear in 7 (see e.g. Fig. 1 in 
Ref. [4]). This is so since in the strict hard core limit 
one expects particles driven with different 7, to follow 
the same path through phase space, only with different 
velocities cx 7, and it then follows that many quan- 
tities (e.g. the forces) are just proportional to 7[4, 5]. 
The alternative is a recently deviced method to perform 
shearing simulations with hard particles [3, 6]. 

The transition in shear-driven systems still appears to 
be rather poorly understood. There is e.g. no accepted 
value for the exponent for the divergence of the viscosity; 
determined values range between 2.0 and 2.8[2, 5, 7-11], 
and this appears to, to at least some extent, be because of 
the lack of understanding of the mechanism behind this 
divergence. To illustrate the complications we point out 
that one typically expects both shear viscosity 77 = cr/7 
and the pressure-equivalent quantity, rjp = p/7, to di- 
verge in the same way, but since /i = a/p = if /rjp has 



a pronounced density dependence [3, 11] in the relevant 
density interval, naive fits of ?]((/') and rjp{(j)) to algebraic 



divergences. 



, give differing values for the criti- 



cal parameters. One way to resolve this issue is to include 
corrections to scaling in the analyses, but even though 
such a program has been successfully accomplished [2], 
this requires very high precision data very close to the 
transition and the scaling analysis becomes both difficult 
and opaque. 

In this Letter we present results from relaxation simu- 
lations which are done by first driving at a certain shear 
rate and then stopping the shearing and letting the sys- 
tem relax according to its dynamics. The relaxation time 
Troiax is then determined from the decay of the energy. 
We believe that this relaxation time is a fundamental 
quantity which is at the root of the divergence of pressure 
and shear viscosity. We also consider another time, Tdiss, 
which is related to the rate at which energy is dissipated 
in steady shearing and show that this quantity behaves 
similarly to Tioiax- We further show that 77 and rjp may 
be written as products of Tdiss and some 0-dependent 
correction factors, and we show that this picture nicely 
explains the need for corrections to scaling in the scaling 
analysis of Ref. [2] . The methods suggested here should 
be useful for studies of the jamming transition through 
both simulations and experiments. 

We simulate frictionless soft disks in two dimensions 
using a bi-dispersive mixture with equal numbers of disks 
with two different radii of ratio 1.4. Length is measured 
in units of the diameter of the small particles {dg = 1). 
With rij the distance between the centers of two parti- 
cles, dij the sum of their radii, and the relative overlap 
5ij = 1 — rij /dij for rij < dij and Sij = otherwise, the 
interaction between the particles is 



e6l/2, 



with e = 1. We use Lees-Edwards boundary conditions 
[12] to introduce a time-dependent shear strain ^(t) ~ ^7. 
With periodic boundary conditions on the coordinates Xi 
and iji in an L x L system, the position of particle i in 
a box with strain 7 is defined as = {xi + jyi, yi). We 
simulate overdamped dynamics at zero temperature with 
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FIG. 1. Relaxation of the energy at difTerent 4>. The fig- 
ure shows the relaxation of energy after the shearing has 
been switched off. The preceeding shearing simulations were 
performed at very low shear rates in order to stay close to 
the linear region; the densities and the initial shear rates 
were {4>,-y) = (0.834,10"*), (0.838,10"*), (0.840,5 x 10"^), 
(0.8408,2 X lO"''), (0.8416,10""). To determine the relax- 
ation times, Treiax, we fit the energy to an exponential decay, 
only using data with E < 10"^^. 



the equation of motion [13], 

with kd — I. In this model dissipation occurs when the 
particles move relative to the steady shearing velocity 
yijx. The effects of instead letting the dissipation be 
given by the relative velocity of particles in contact will 
be discussed elsewhere[14]. 

The key quantity in this Letter, the relaxation time, is 
determined through a two-step process: first the system 
is driven in steady shear at a constant shear rate 7; then 
the shearing is stopped and the system is allowed to relax 
down to a minimum energy. As the simulations discussed 
here are at densities somewhat below (jjj, the final state is 
always a state of zero energy, and after a short transient 
time, the decay is exponential, 

E{t) exp(-t/Ti.clax)- 

A few realisations of such relaxations are shown in Fig. 1. 
For each realisation the relaxation time is determined 
from the data with E{t) < 10"^^, where the decay is ex- 
ponential to an excellent approximation. We determine 
TVciaxl^j 7) as the average relaxation time from about 10- 
100 such relaxations. 

Figure 2(a), which is Trciax versus </> in a narrow den- 
sity interval just below 0j, shows that Trdax increases 
very rapidly with 0. The data are shown for a few differ- 
ent 7 and we conclude that Trdax at a given approaches 
a well-defined limiting value as 7 decreases and the linear 
region is approached. In this Letter we analyze the data 
within (or close to) this linear region only; the behavior at 
larger 7 will be examined elsewhere. We first determine 
the critical behavior from the eight points in Fig. 2(a) 
which are in the linear region and close below jamming, 
i.e. the points with the lowest shear rate for each density 
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FIG. 2. Behavior of the relaxation time. Panel (a) shows how 
TVciax depends on both <f> and 7, which is here the shear rate 
of the preparatory run (the relaxations are performed with 
7 = 0). Each value here is the average of relaxation times de- 
termined from of a large number of different relaxations. At 
sufficiently low 7, Treiax approaches well-defined values that 
only depend on (j>. Panel (b) is a double-log plot, only includ- 
ing the points with small enough 7 to be in the linear region. 
The figure shows a fit to the eight points with (j> > 0.834; the 
points with (j)j — (j> > 0.01 are not included in the fit. 



in the range 0.834 < <j) < 0.8416. Fitting these points 
to an algebraic divergence, Ti.ciax(0) = ^|<50|^'', (where 
S(j> = (j>j — (/)) gives Fig. 2(b) and the critical parameters 
= 0.8433 ± 0.0002 and /3 = 2.73 ± 0.15, which arc in 
good agreement with Refs. [2, 15]. The quoted errors rep- 
resent max/min-values, corresponding to three standard 
deviations in the estimated quantities. 

Our assumption is that this increase of the relaxation 
time as jamming is approached is the fundamental phe- 
nomenon which is at the root of the divergence of other 
quantities as e.g. the shear viscosity, rj. The relaxation 
mode should be related to the isolated mode with fre- 
quency Wmin in Ref. [3], and we expect Treiax ^m\n- 

(The different powers of time here reflect the differences 
in dynamics. In overdamped dynamics there is a velocity 
that is proporional to a force, whereas one in vibrational 
analyses assumes Newtonian dynamics with massive par- 
ticles where the acceleration is proportional to the force.) 
Note also that the lowest mode being isolated explains 
why the relaxation is almost perfectly exponential after 
the initial decay. One would otherwise typically expect 
E{t) to be given by a sum of several modes with close but 
different time constants. A further result from Ref. [3, 16] 
is that (and thereby Trdax) diverges with the same 

exponent as the shear viscosity. This conclusion will also 
be reached below in a different way. 

We will now try and establish a link between the shear 
viscosity, rj, which is typically measured in both simu- 
lations and experiments, and the above obtained Treiax- 
This will be done in two steps: we first derive an expres- 
sion for a similar time, Tdissj which is obtained directly in 
the shear driven simulations; we then express rj in terms 

of Tdiss- 

The expression for the dissipation time, Tdiss, is ob- 
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FIG. 3. Dissipation time Tdiss obtained in shear driven simu- 
lations. Panel (a) is rdiss against (j) at different 7. Panel (b) 
shows the data considered to be in the linear region against 
4>.j — 4>- The dashed line shows the result of a fit to the eight 
points with (pj - (p < 0.01, giving (pj = 0.8434 and /3 = 2.84. 



tained from a power balance. The idea is that the sup- 
pUed power, which is 177 per unit area, on average should 
be balanced by the dissipated power. Defining Tdiss such 
that E/ Tdiss is the rate at which the energy is dissipated 
defines 



E 

(T7 



(1) 



Note that it follows directly that Tdiss should diverge 
with the exponent /3 since a/j ^ \S(j)\~^ and E/j^ ^ 
|<5<A|-2/3[4].^ 

The dissipation time Tdiss versus (j) for a few different 
shear rates is shown in Fig. 3(a). Just as for Tj-oiax we find 
that Tdiss increases rapidly when (f) increases towards (t>j, 
and we also find well-defined low-7 limits, with deviations 
for larger 7. Here Tdiss becomes smaller for larger 7, 
which is the same behavior as in the shear viscosity, but 
opposite to the behavior of Trdax, discussed above. 

The rational to introduce Tdiss was to find a quantity in 
the shearing simulations that behaves similarly to Ti-ciax, 
and thus establish a link between the relaxation dynam- 
ics and the shearing simulations. It is however clear that 
these two quantities cannot be identical. Since the ini- 
tial dissipation in a relaxation simulation has to be the 
same as the dissipation under steady shear, Tdiss is equal 
to the initial decay rate in a relaxation simulation. Tj-ciax 
on the other hand is the decay rate at long times. This 
means that Tdiss should get contributions from all de- 
cay modes that are present in the system. Ti-dax, on the 
other hand, is determined by the slowest mode only, since 
that is the only mode that persists after sufficiently long 
times. Since Tdiss gets contributions from modes with 
smaller time constants it follows that Tdiss < Troiax- This 
is confirmed by Fig. 4 which, furthermore, shows that 
"Jdiss /'''relax iucrcascs with increasing (/> and appears to 
approch unity as — >■ </)j. We relate this to the ob- 
servation in Ref. [3] that the relative contribution to the 
shear viscosity of the isolated mode (in their notation, 
(To/ct) approaches unity as jamming is approached, which 
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FIG. 4. Relation between dissipation time and relaxation 
time. The figure shows that Tdiss /Trdax increases slowly as 
(pj is approached from below. The solid line shows a fit to 
an algebraic function which approaches unity at ^j. The 
dashed line is a parametrization to help compare the size of 
this correction with the other corrections to scaling in Fig. 5. 



means that the weight of the other modes decreases. We 
likewise expect the contributions from the faster modes 
to Tdiss to become less important as 0j is approached, 
which implies Tdiss /TVeiax 1- We summarize the above 
in terms of two conclusions of importance for the present 
work: (i) Properties determined in steady shear will nec- 
essarily be different from the properties determined from 
the long-time behavior of the relaxation simulations, (ii) 
This difference is however rather small and one should 
therefore expect results based on T^oiax and Tdiss, respec- 
tively, to be very similar. 

A fit of Tdiss to the algebraic divergence (where we 
again use only the eight points in the linear region and 
close to is shown in Fig. 3(b) and gives = 0.8434± 
0.0003 and (i = 2.84±0.20. Both values are close to (just 
slightly higher than) the corresponding values from the 
analysis of Trciax, and this again suggests that Tdiss is a 
good approximation of Trdax ■ 

The finding that Tdiss to a good approximation diverges 
algebraically, gives a ground for understanding the need 
for corrections to scaling in the analyses of t] and in 
Ref. [2]. For the 7—^0 limit at densites below </)j, cor- 
rections to scaling means that the divergence cannot be 
well approximated by the algebraic ^jiJ^i'l^'^ alone, but 
that one instead has to use 



A|,50r^ (1 + ai5(/.r 



(2) 



which follows from using h — in the unnumbered 

equation before Eq. (3) in Ref. [2]. Here the correction 
to scaling exponent to appears together with the corre- 
lation length exponent v. This behavior is illustrated in 
Fig. 5(a) which shows both 77 and t]^ versus 0j — (/). Also 
shown is r\E — Ve/'^ which behaves the same as rjp, to 
an excellent approximation. (This is so since p ~ '^ij ^ij 
whereas E ~ ^j^- (5fj [17].) As is clear from the figure, 
and rjp behave differently, and attempts to determine /3 
from algebraic fits without corrections, give P = 2.35 and 
2.59, respectively, as shown by the solid lines. Since one 
expects the asymptotic behavior of ry and rjp to be the 
same, this discrepancy calls for including corrections to 
scaling (as in Eq. (2)), which was also done successfully 
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FIG. 5. Viscosities r], r/p, and tie in the light of Eqs. (3) 
and (4). The dashed lines in panel (a) show the divergence 
of Tdiss with /? = 2.84. The symbols show how 77 and rjp 
approach the presumed asymptotic scaling behavior, and it 
is clear that the corrections to scaling are much larger for rj 
than for rip. The same is seen by naively fitting 77 and r]p to 
algebraic divergences (shown by solid lines) which give /3 — 
2.35 and 2.59, respectively, where we note that the exponent 
obtained for rj is further off the asymptotic value P = 2.84. 
Panel (b) shows both fiE ~ ct/^/E, the similar dimensionless 
friction /i = a/p, and (the correction in Eq. (3)). This last 
quantity appears to be linear in (f), which is consistent with 
uji/ ~ 1 as found in Ref. [2] . 



in Ref. [2]. 

We will now shovir that 77 and 77^ may be written as 
products of Tdiss and some correction factors. After in- 
troducing He ~ o I \fE in analogy with the dimensionless 
friction fx — a/p (see Fig. 5(b)) we find using Eq. (1) that 



riE = VE/j = He Tdis 



(3) 
(4) 



We note two things: (i) That the corrections of 77 and 
TjE are and he, respectively gives a very direct ex- 
planation to why the correction to scaling in rjp (which 
behaves essentially the same as tie) is so much smaller 
than in 77(2]. See also below for a direct comparison of 
these correction terms, (ii) As shown in Fig. 5(b) h^ (the 
correction factor that goes together with rj) is linear in 
(j) to an excellent approximation and the same holds for 
He-, though in a more narrow range below (/),;. Together 
with Eq. (2) we therefore conclude that luv w 1, again in 
agreement with Ref. [2]. 

As discussed above our starting assumption is that 
Treiax diverges algebraically and it then follows from Fig. 4 
that Tdiss is given by this algebraic divergence times a cor- 
rection factor. To argue that the critical behavior should 
be determined from Tdiss rather than 77 or 77^, we now 
want to show that this correction in Tdiss that one can- 
not eliminate (if one only has access to data from steady 
shearing) is considerably smaller than the correction fac- 
tors in 77 and 77^;. To do that we write each correction on 
the form (l+a|(50|) and compare the magnitude of "a" for 



the different cases. We then find h% ~ O.OO35(H-228|50|) 
and (close to 0j) he ~ 0.061(1 -I- 88|50|) and note that 
both these correction terms are clearly bigger than the 
correction in Tdiss/Trciax ~ (1 — 17|(5(/)|) from Fig. 4[18]. 
This strengthens our confidence in the use of Tdiss for 
determining the critical behavior, though it is of course 
Trciax that is the ideal quantity for such analyses. 

The results above should also be useful for analyzing 
experiments, but instead of using Tdiss — E/a-j one could 
then make use of Tdiss ^ /crj which is an expression in 
terms of pressure instead of the elastic energy. This could 
be advantageous since pressure should be more readily 
available in experiments than energy. 

The relaxation dynamics around the jamming tran- 
sition has been studied before, but then with a rather 
different preparation of the starting configurations [17]. 
In that study configurations were first generated ran- 
domly, then relaxed to a zero-energy state with the con- 
jugate gradient method, and after that perturbed by a 
pure affine shear deformation. The relaxation time was 
then determined from the relaxation of such initial states 



by fitting the shear stress to cr{<p,t) ^ t "e 



-t/7 



with 



a — 0.55(5), and was found to diverge as t ^ {<j)j — 
with C = 3.3(1). This exponent is clearly bigger than 
our /3 = 2.73 ± 0.15. One possible explanation for this 
difference is that in the present study we have been very 
careful to apply a slow shear driving in the preparation 
step, whereas they in their work apply the pure shear 
deformation suddenly, which should be more like a rapid 
shearing. Indeed, as shown in Fig. 2(b) any given fixed 
shear rate would give too large values for Tj-eiax as one gets 
close to (f)j, and from analyses of such data one would 
expect to get a too high value of the exponent for the 
divergence. 

To conclude, we have determined Treiax from relax- 
ational simulations and suggest that the slowing down 
of the relaxation as is approached is the fundamental 
reason for the divergence of rj and other similar quanti- 
ties. Strong support for this idea is obtained from the 
finding by others that there is an isolated mode that 
dominates the behaviour close to (/)j[3]. We have fur- 
ther introduced Tdiss which is determined directly in shear 
driven simulations and have shown that these two quanti- 
ties, in the linear region and close to are very similar. 
From the connection between Tdiss and rj we further argue 
that the need for corrections to scaling in analyses of 77 
and related quanties is largely due to the ^-dependence 
of He = (j/VE. Our results should also be helpful for 
getting more accurate determinations of the critical be- 
haviour from experimental data. 
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